﻿
#include <iostream>
#include <math.h>

static const double pi = (3.14159265358979323846);

/// calculates 1*2*3*...*n
inline double factorial(int n);

/// calculates 1*3*5*...*2n-3
inline double factorial2(int n);

/// calculates int(cos(at)dt)^n from t1 to t2;
inline double cosIntN(int a, double t1, double t2, int n);

/// Calculates the arclength of elliptic curve starts from theta1 and ends at 
/// theta2, using infinite series approach given in the reference below 
/// NOTE:
/// if theta1 > theta2 arc length will be negative
/// maximum number of iterations is set as 100, change from code if necessary
/// Inputs:
/// theta1 and theta2 are defined such that;
/// defined in ellipse parameter such that for a point (x,y) on ellipse:
/// x = a*cos(theta1)
/// y = b*sin(theta2)
/// a,b are semimajor and semiminor axes of ellipse, however function
/// supports entering them as vice versa
/// tolerance sets the accuracy of calculation
/// Outputs:
/// arcLen: arclength in units of the input a and b
/// precision: difference of ellipse length with one higher iteration
/// calculation
/// n: number of iterations
/// referenced from:
/// http://pages.pacificcoast.net/~cazelais/250a/ellipse-length.pdf
/// example:
/// [arcLen, precision, n] = ellipseArcLength(0, 2*pi, 5,4,1e-5)
/// Author: Mehmet Burak Ekinci
/// Mail: elessar208@gmail.com
inline double ellipseArcLength(double theta1, double theta2, double a, double b, double tolerance,
    double &precision, int &n)
{
    double arcLen;

    if (b > a) {  // compensate for elliptic integral a>b constraint
      theta1 = theta1 - pi / 2;
      theta2 = theta2 - pi / 2;
      double temp = a;
      a = b;
      b = temp;
    }

    n = -1;
    arcLen = 2 * (theta2 - theta1);
    double eps2 = -(1 - (b*b) / (a*a));
    double c;
    //  epsilon squared 
    double q = tolerance + 1;
    while (abs(q) > tolerance ) {
      n = n + 1;  // n will start from 0 this way
      c = cosIntN(1, theta1, theta2, 2 * n);
      q = pow((-1) , (n + 1)) * factorial2(n) / pow(2,  n) / factorial(n) * pow(eps2, n) * c;
      arcLen = arcLen + q;
      if (n > 100) {
        std::cout << "series did not converge in 100 iterations" << std::endl;
        break;
      }
    }
    arcLen = arcLen * a;
    precision = q;

    return arcLen;
}

/// calculates 1*2*3*...*n
inline double factorial(int n)
{
    double f = 1;
    if (n <= 1) {
        f = 1;
        return f;
    }
    for (int i = 1; i <=n; i++) {
       f = f * i;
    }
    return f;
}

/// calculates 1*3*5*...*2n-3
inline double factorial2(int n)
{
    double f = 1;
    if (n <= 1) {
        f = 1;
        return f;
    }
    for (int i = 2; i <=n; i++) {
       f = f * (2 * i - 3);
    }
    return f;
}

/// calculates int(cos(at)dt)^n from t1 to t2;
inline double cosIntN(int a, double t1, double t2, int n)
{
    double I = 0;
    if (n == 0) {
        I = t2 - t1;
        return I;
    }

    double A = pow(cos(a*t2), (n-1))*sin(a*t2)/n/a - pow(cos(a*t1),(n-1))*sin(a*t1)/n/a;
    if (n == 1) {
        I = A;
    } else {
        double B = double(n-1)/double(n)*cosIntN(a,t1,t2,n-2);
        I = A + B;
    }
    return I;
}

int main(int argc, char const *argv[])
{

    std::cout.precision(15);
    std::cout.unsetf(std::ios::scientific); // 取消科学计数法
    std::cout.setf(std::ios::fixed); //强制显示小数点后的无效0


    double precision = 0;
    int n = 0;

    // double I = cosIntN(1, 0.0, 2.0*pi, 2);
    // std::cout << "cosIntN:" << I << std::endl;

    double arcLen = ellipseArcLength(0.0, 2.0*pi, 10, 6, 1e-10, precision, n);

    std::cout << "arcLen:" << arcLen << "  precision: " << precision << "  n:" << n << std::endl;

    double totalLen = 0.0;
    for (int i=0; i< 27; i++) {
       arcLen = ellipseArcLength(double(i)*2*pi/27, double(i+1)*2*pi/27, 10, 6, 1e-10, precision, n); 
       // std::cout << "no:" << i << " rad:" << double(i)*2*pi/27 << " arcLen:" << arcLen <<std::endl ;
       std::cout << "no:" << i << " arcLen:" << arcLen <<std::endl ;
       totalLen += arcLen;
    }
     std::cout << "sum all:" << totalLen <<std::endl ;

    return 0;
}
